home *** CD-ROM | disk | FTP | other *** search
/ CD Ware Multimedia 1995 May / cd Ware (Juegos) Epimundo.iso / DOS / C / CMATH.ZIP / LOG10.C < prev    next >
Encoding:
C/C++ Source or Header  |  1989-03-21  |  4.8 KB  |  227 lines
  1. /*                            log10.c
  2.  *
  3.  *    Common logarithm
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, log10();
  10.  *
  11.  * y = log10( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns logarithm to the base 10 of x.
  18.  *
  19.  * The argument is separated into its exponent and fractional
  20.  * parts.  The logarithm of the fraction is approximated by
  21.  *
  22.  *     log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
  23.  *
  24.  *
  25.  *
  26.  * ACCURACY:
  27.  *
  28.  *                      Relative error:
  29.  * arithmetic   domain     # trials      peak         rms
  30.  *    IEEE      0.5, 2.0     30000      1.5e-16     5.0e-17
  31.  *    IEEE      0, MAXNUM    30000      1.4e-16     4.8e-17
  32.  *    DEC       1, MAXNUM    14000      2.5e-17     6.0e-18
  33.  *
  34.  * In the tests over the interval [1, MAXNUM], the logarithms
  35.  * of the random arguments were uniformly distributed over
  36.  * [0, MAXLOG].
  37.  *
  38.  * ERROR MESSAGES:
  39.  *
  40.  * log10 singularity:  x = 0; returns MINLOG
  41.  * log10 domain:       x < 0; returns MINLOG
  42.  */
  43.  
  44. /*
  45. Cephes Math Library Release 2.1:  December, 1988
  46. Copyright 1984, 1987, 1988 by Stephen L. Moshier
  47. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  48. */
  49.  
  50. #include "mconf.h"
  51. static char fname[] = {"log10"};
  52.  
  53. /* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
  54.  * 1/sqrt(2) <= x < sqrt(2)
  55.  */
  56. #ifdef UNK
  57. static double P[] = {
  58.   4.58482948458143443514E-5,
  59.   4.98531067254050724270E-1,
  60.   6.56312093769992875930E0,
  61.   2.97877425097986925891E1,
  62.   6.06127134467767258030E1,
  63.   5.67349287391754285487E1,
  64.   1.98892446572874072159E1
  65. };
  66. static double Q[] = {
  67. /* 1.00000000000000000000E0, */
  68.   1.50314182634250003249E1,
  69.   8.27410449222435217021E1,
  70.   2.20664384982121929218E2,
  71.   3.07254189979530058263E2,
  72.   2.14955586696422947765E2,
  73.   5.96677339718622216300E1
  74. };
  75. #endif
  76.  
  77. #ifdef DEC
  78. static short P[] = {
  79. 0034500,0046473,0051374,0135174,
  80. 0037777,0037566,0145712,0150321,
  81. 0040722,0002426,0031543,0123107,
  82. 0041356,0046513,0170752,0004346,
  83. 0041562,0071553,0023536,0163343,
  84. 0041542,0170221,0024316,0114216,
  85. 0041237,0016454,0046611,0104602
  86. };
  87. static short Q[] = {
  88. /*0040200,0000000,0000000,0000000,*/
  89. 0041160,0100260,0067736,0102424,
  90. 0041645,0075552,0036563,0147072,
  91. 0042134,0125025,0021132,0025320,
  92. 0042231,0120211,0046030,0103271,
  93. 0042126,0172241,0052151,0120426,
  94. 0041556,0125702,0072116,0047103
  95. };
  96. #endif
  97.  
  98. #ifdef IBMPC
  99. static short P[] = {
  100. 0x974f,0x6a5f,0x09a7,0x3f08,
  101. 0x5a1a,0xd979,0xe7ee,0x3fdf,
  102. 0x74c9,0xc66c,0x40a2,0x401a,
  103. 0x411d,0x7e3d,0xc9a9,0x403d,
  104. 0xdcdc,0x64eb,0x4e6d,0x404e,
  105. 0xd312,0x2519,0x5e12,0x404c,
  106. 0x3130,0x89b1,0xe3a5,0x4033
  107. };
  108. static short Q[] = {
  109. /*0x0000,0x0000,0x0000,0x3ff0,*/
  110. 0xd0a2,0x0dfb,0x1016,0x402e,
  111. 0x79c7,0x47ae,0xaf6d,0x4054,
  112. 0x455a,0xa44b,0x9542,0x406b,
  113. 0x10d7,0x2983,0x3411,0x4073,
  114. 0x3423,0x2a8d,0xde94,0x406a,
  115. 0xc9c8,0x4e89,0xd578,0x404d
  116. };
  117. #endif
  118.  
  119. #ifdef MIEEE
  120. static short P[] = {
  121. 0x3f08,0x09a7,0x6a5f,0x974f,
  122. 0x3fdf,0xe7ee,0xd979,0x5a1a,
  123. 0x401a,0x40a2,0xc66c,0x74c9,
  124. 0x403d,0xc9a9,0x7e3d,0x411d,
  125. 0x404e,0x4e6d,0x64eb,0xdcdc,
  126. 0x404c,0x5e12,0x2519,0xd312,
  127. 0x4033,0xe3a5,0x89b1,0x3130
  128. };
  129. static short Q[] = {
  130. 0x402e,0x1016,0x0dfb,0xd0a2,
  131. 0x4054,0xaf6d,0x47ae,0x79c7,
  132. 0x406b,0x9542,0xa44b,0x455a,
  133. 0x4073,0x3411,0x2983,0x10d7,
  134. 0x406a,0xde94,0x2a8d,0x3423,
  135. 0x404d,0xd578,0x4e89,0xc9c8
  136. };
  137. #endif
  138.  
  139. #define SQRTH 0.70710678118654752440
  140. #define L102A 3.0078125E-1
  141. #define L102B 2.48745663981195213739E-4
  142. #define L10EA 4.3359375E-1
  143. #define L10EB 7.00731903251827651129E-4
  144.  
  145. extern double LOGE2, SQRT2, MAXLOG, MINLOG;
  146.  
  147. double log10(x)
  148. double x;
  149. {
  150. int e;
  151. short *q;
  152. double w, y, z;
  153. double frexp(), ldexp(), polevl(), p1evl();
  154.  
  155. /* Test for domain */
  156. if( x <= 0.0 )
  157.     {
  158.     if( x == 0.0 )
  159.         mtherr( fname, SING );
  160.     else
  161.         mtherr( fname, DOMAIN );
  162.     return( MINLOG );
  163.     }
  164.  
  165. /* separate mantissa from exponent */
  166.  
  167. #ifdef DEC
  168. q = (short *)&x;
  169. e = *q;            /* short containing exponent */
  170. e = ((e >> 7) & 0377) - 0200;    /* the exponent */
  171. *q &= 0177;    /* strip exponent from x */
  172. *q |= 040000;    /* x now between 0.5 and 1 */
  173. #endif
  174.  
  175. #ifdef IBMPC
  176. x = frexp( x, &e );
  177. /*
  178. q = (short *)&x;
  179. q += 3;
  180. e = *q;
  181. e = ((e >> 4) & 0x0fff) - 0x3fe;
  182. *q &= 0x0f;
  183. *q |= 0x3fe0;
  184. */
  185. #endif
  186.  
  187. /* Equivalent C language standard library function: */
  188. #ifdef UNK
  189. x = frexp( x, &e );
  190. #endif
  191.  
  192. #ifdef MIEEE
  193. x = frexp( x, &e );
  194. #endif
  195.  
  196. /* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
  197.  
  198. if( x < SQRTH )
  199.     {
  200.     e -= 1;
  201.     x = ldexp( x, 1 ) - 1.0; /*  2x - 1  */
  202.     }    
  203. else
  204.     {
  205.     x = x - 1.0;
  206.     }
  207.  
  208.  
  209. /* rational form */
  210. z = x*x;
  211. y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
  212. y = y - ldexp( z, -1 );   /*  y - 0.5 * x**2  */
  213.  
  214. /* multiply log of fraction by log10(e)
  215.  * and base 2 exponent by log10(2)
  216.  */
  217. z = (x + y) * L10EB;  /* accumulate terms in order of size */
  218. z += y * L10EA;
  219. z += x * L10EA;
  220. w = e;
  221. z += w * L102B;
  222. z += w * L102A;
  223.  
  224.  
  225. return( z );
  226. }
  227.